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We investigate a system of fermions on a two-dimensional optical square lattice in the strongly 
repulsive coupling regime. In this case, the interactions can be controlled by laser intensity as well 
as by Feshbach resonance. We compare the energetics of states with resonating valence bond d- 
wave superfluidity, antiferromagnetic long range order and a homogeneous state with coexistence 
of superfluidity and antiferromagnetism. We show that the energy density of a hole Choieix) has a 
minimum at doping x = Xc that signals phase separation between the antiferromagnetic and d-wave 
paired superfluid phases. The energy of the phase-separated ground state is however found to be 
very close to that of a homogeneous state with coexisting antiferromagnetic and superfluid orders. 
We explore the dependence of the energy on the interaction strength and on the three-site hopping 
terms and compare with the nearest neighbor hopping t-J model. 

PACS numbers: 37.10.Jk, 74.72.-h, 74.20.-z 



I. INTRODUCTION 



The basic mechanism by which attractively interacting 
fermions can mrdcrgo Bardeen-Coopcr-Schrieffer(BCS)- 
typc transition through s-wave pairing is now well under- 
stood. The copper-oxide (CUO2) materials are the first 
example of superconductivity arising from strong repul- 
sive interactions It is argued that the single-band 
Hubbard model captures the essential physics of on-site 
repulsion. From variational calculations of the Hubbard 
model, the mechanism is now understood to be an an- 
tiferromagnetic exchange mechanism that favors singlet 
pairs on different sites leading to pairing in the d-wave 
channel 0- These calculations show a strong deviation 
from the standard BCS paradigm with a separation of 
energy scales for pairing and long range coherence 0, |3] . 
However, in the absence of a rigorous solution of the Hub- 
bard model in two dimensions, the nature of the ground 
state, and in particular, whether it is a d-wave supercon- 
ductor or not is still in debate Q. 



Given the difhculties of solving the Hubbard model, an- 
other route that is being currently attempted is to turn 
to ultracold atoms in optical lattices 0, 0] and use them 
to emulate Hubbard-like models long studied in the con- 
densed matter physics. Unlike the materials studied in 
the condensed matter physics in which the interaction 
strengths are usually fixed, in optical lattices the inter- 



actions can be varied over a wide range from weak to 
strong coupling. Furthermore, in the condensed matter 
systems changing carrier concentration usually modifies 
the degree of disorder, while the optical lattices are per- 
fect and clean. Thus, the optical lattices allow direct 
simulation of one quantum system by another Q . It is 
hoped that a deeper understanding of superconductivity 
in repulsive models will suggest mechanisms for driving 
the transition temperature higher. At present the main 
bottlenecks faced in creating low entropy states in opti- 
cal lattices pertain to cooling mechanisms and the role of 
inhomogcneity introduced by an overall trapping poten- 
tial. Also, the multi-band effects when strongly interact- 
ing fermions are loaded in the optical lattice have to be 
carefully taken into account 



le opt K 



In this article we discuss the problem of repulsively in- 
teracting fermions in a square optical lattice. The repul- 
sion can be varied either by increasing the laser intensity 
or by a Feshbach resonance. For both cases we derive 
an effective Hamiltonian in the strong coupling regime 
that is independent of the interaction mechanism. We ex- 
pect both d-wave superconducting and antiferromagnetic 
fluctuations to dominate in the ground state. The main 
question we ask is whether the ground state energetically 
chooses to phase separate or to stay in a homogeneous 
state as doping and interactions are changed. In the 
past, analogous questions were addressed in perovskite 
compounds such as manganites In that case, 

the phase separation between undoped antiferromagnetic 
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and hole-doped ferromagnetic regions was observed. The 
origin of the interaction was, however, coulombic repul- 
sion between the localized and mobile electrons. 



Our result is that the Hubbard model does show a ten- 
dency to phase separate at low hole density. This result is 
based on a variational theory that compares the energet- 
ics for two possible ground states: one that allows phase 
separation and another with homogeneous coexistence. 
Here, the possible biases are applied on equal degree to 
these different scenarios. Given the small difference in en- 
ergy between the two solutions for phase separation and 
homogeneous coexistence, other factors may determine 
the ultimate behavior of the system: for example, long- 
range interactions may tend to favor coexistence whereas 
disorder may drive the system toward phase separation. 
Also, the effect of trapping potential needs to be better 
understood [l3.[l5}. 



II. MODEL 



Hamiltonian 



The interactions between fermions in the optical lat- 
tices can be controlled in two ways 0, : 



i) Laser Intensity controlled: By varying the laser in- 
tensity or the optical lattice depth Vo, we can mod- 
ify the tunneling rate t of atoms between neighbor- 
ing sites and the on-site interaction U. t and U can 
be calculated from the band structure theory. In the 
deep lattice limit, the tunneling rate is given hy t k. 



3/4 



action U 



exp 



—akpEr 



1/2 



3/4 



and the on-site inter- 



where a is the s- 



wave scattering length that is fixed and kp is the Fermi 
wave vector. For realizing repulsive fermions, potassium 
isotope "^^K can be used where a for hyperfinc states is 
positive (« 160ao [HI) with defined as Bohr radius. 
Typical a for ^Li is w — 2000ao. Here = 2mA^ 
recoil energy of the optical lattice generated by a set of 
counter-propagating laser beams of wavelength A. The 
expressions above are valid in the non-resonant regime 
where kpa 1. It is then possible to restrict the Hilbert 
space to a single band Hubbard model given by 

<i,j><J i 

where a — ±1 is the pseudo-spin index. The sums 



run over all pairs of nearest-neighbor lattice sites and 
TT-i^cr = c\^Ci^rT- Wc define doping parameter a; = 1 — 
N/Niatt where is the total number of fermions of both 
spin species and Niatt is the number of sites in the optical 
lattice. 



ii) Feshbach resonance controlled : By tuning the scat- 
tering length a using a magnetic field B. Scatter- 
ing length near resonance takes the form a = a6g[l ~ 
AB/{B — Bo)] where abg is the background scatter- 
ing length (« leOflo for ^°K). Also AB « 8G and 
Bo « 202 G for "^^K. Here, the scattering length a can 
be tuned from the positive values (repulsive) to zero (non- 
interacting) and to the negative values(attractive). Near 
the resonance, the interactions can be easily driven into 
a regime where they exceed the band gap or the energy 
difference between the lowest two levels in a single lat- 
tice site. In this resonant limit the wave function of the 
two-particle state includes contributions from all excited 
states and does not reduce to a symmetrized product 
of the two single particle ground state wave functions. 
Thus, the minimal set of basis per site should include a 
single particle ground state as well as a state with dou- 
ble occupancy. Consequently, the overlap of the wave 
functions for nearest neighbor sites and the hopping am- 
plitudes depend on the particular configuration. This 
situation is illustrated in the Fig. [TJ The physics of 
the resonant two-particle states introduces two additional 
parameters St and Sg, which are incorporated into the 
Hamiltonian in the following form 



mcorpor 

Mm 



<i,j>a 

+Uo^rii,+ini^_i. (2) 



Here, t is controlled by the optical lattice depth Vo as in 
Eq. [T] However in the Feshbach resonance case, not just 
the optical lattice depth Vo but also applied magnetic 
field and other physical parameters of the system deter- 
mine the value of Uo- Further discussions below and in 
the Appendix show that within second order perturba- 
tivc theory, the Hamiltonians of Eq. [l]and Eq. [2] become 
analogous with the identification of an effective coupling 
U. 



B. Our Approach 



Given the Hamiltonian, the next step is to determine 
the quantum phases arising from the strongly repulsive 
interactions among fermions in two dimensions. At half 
filling, or zero doping {x = 0), the system is a Mott 
insulator with a finite gap in the charge spectrum and 
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FIG. 1: (Color online) Here, we illustrate the effect of having a resonant two-particle state away from the single particle 
ground state. The hopping amplitudes are essentially given as the overlap of the wave functions of the neighboring sites. Thus, 
within this "two level" model we can identify three different hopping processes and amplitudes, a) from a single particle state 
to another single particle state, b) from a two-particle state to a single particle state or the reverse process, and c) from a 
two-particle state to another two-particle state. Blue and red colors indicate fermions with opposite spins. 



has long range antiferromagnetic order. As the system is 
doped with holes, antiferromagnetic order decreases and 
the system develops superfluidity described by a coher- 
ent superposition of the resonating valence bond (RVB) 
singlets Here, each singlet is formed by a pair of 
opposite fermion species in an orbital state with d-wave 
spatial symmetry. The question that arises is: What are 
the characteristics of the quantum phase at small hole 
density x? 

In this article, we discuss two possible scenarios: 

1) Phase Separation: In this scenario, the system is 
macroscopically phase separated into regions of antifer- 
romagnetism and superfluidity. We calculate the behav- 
ior of the energy density of a hole ehoie{x) to check for 
phase separation. 

2) Homogeneous: Here, we consider a single homoge- 
neous phase with coexisting antiferromagnetic and super- 
fluid orders occupying the whole system. We evaluate the 
energetics of variational states that include both super- 
fluid and antiferromagnetic orders. 

By comparing with the more restricted standard t-J 
model we find that the presence of effective 3-site hopping 
terms in the full Hamiltonian defined in Eq. |4] leads to a 
lower energy at non-zero hole density x. We also find a 



stronger tendency to phase separate in the full Hubbard 
Hamiltonian compared to the more restrictive t-J model. 



C. Effective Hamiltonian in Strong Coupling 



In the large U/t limit, the Hamiltonian in Eq. [T]can 
be transformed into a block diagonal form where the 
blocks preserve the total number of the doubly occupied 
sites and the different blocks are connected by the ki- 
netic/hopping operators (see Appendix). The effective 
Hamiltonian T-L^f f in the stron g co upling limit is obtained 
by a canonical transformation|2ll |22| 

He// = V + /Co + ^[/Ci,/C_i]+0(l/C/2) (3) 

where the kinetic operator ICm changes the number of 
doubly occupied sites by an integer number m. V is 
the interaction term. The term given by a commutator 
(second order term) in Eq. [3] contains two consecutive 
hoppings. Since we consider U/t >> 0, we impose the 
non-double occupancy constraint at each site in order to 
obtain the low energy Hamiltonian: Heff — t- Vc'Hef fVc 
where Vq = 0(1 ~ '^i, 1*^2, 2) is the full Gutzwiller projec- 

i 

tor. As a result, the interaction term V is projected out 
and the Hamiltonian can be cast as 



j.2 

<i,j>cr 



u 



(4) 



<?,z>(T;<2,^>cr 



The summations run over all independent indices and to distinguish it from the more restricted t-J model. In 
Eq. |4] is what we call the projected or full Hamiltonian 
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the second sum, we can identify two-site (when i' = j ) 
and three-site (when i' ^ j) processes. In the i- J model, 
an additional restriction i' = j is imposed such that only 
two-site exchange is allowed. In Figs. [5] and [71 we show 
the effects of this constraint. We have introduced defi- 
nitions: fj^ = hi^ac\ a- ftic creation operator of a singly 

occupied site with spin a particle and d\ ^ = c] ^hi^a-ni^a 
the creation operator of a doubly occupied site by 
adding a particle of spin cr to a singly occupied site with 
a spin a particle where ht^a = 1 — rii^^. Rewriting the 
Hamiltonian in terms of the pseudo-particle operators 
/i,^, fj„, di,a and d\ ^ makes clear the permitted 
physical processes. However, these operators do not 
satisfy the usual anti-commutator nor commutator 
relations and straightforward theoretical treatment is 
difhcult. Thus, we rely on the variational formalism 
whose accuracy depends on the chosen ansatz [23l - [25j . In 
the following section, wc discuss the different variational 
wave functions. 



For the fcrmions near resonance (Eq. [2]) . the same pcr- 
turbative Hamiltonian (Eq. ^ is derived with an effective 
interaction U which has two contributions 



U ■■ 



Uo 



{l + 5g/t) 



2 ■ 



(5) 



Here, Uo is controlled by the optical lattice depth and 
the magnetic field while 5g/t parameter accounts for the 
appearance of a two-particle state that goes beyond the 
single band model (see Appendix). The 5g dependence 
comes from the allowed virtual processes while the 5t 
contribution does not appear in the second order pertur- 
bation. As a consequence, we now consider the phase 
diagram of the Hamiltonian in Eq. |3]as a function of U /t 
and doping x, regardless of the different ways of tuning 
the interaction U . 



D. Variational Wave functions 



1) Phase separated solution using Maxwell con- 
struction: Phase separation results from thermody- 
namic considerations that involve two pure phases. In 
the following, we give their description leaving the details 
of Maxwell construction to the section IIIII (Discussion of 
Results) . 



The pure superfluid variational state is described by 



JVg 



k 



N/2 



|0) 



(6) 



\^ spf) is also called the RVB wave function when Jas- 
trow correlation is set to unity. But as we discuss be- 
low, Y^sFd) is more general and can also describe a 
system with long range antiferromagnetism for an ap- 
propriate choice of the Jastrow correlations. \SFd)MF 
is the solution of the simple mean field Hamiltonian of 
the form: Umf = Eka ^kc[ ,,Ck,a -I- Ek ^k[c[ .^cLk,; + 
h.c.]. Thus, the parameters ak arc given as ak = 

Ak/ 



The d-wave parametrization 



Ak = Arf(cos fc^ — cos ky) is always energetically favorable 
in comparison with the s-wave symmetric parametriza- 
tion Ak = A^. We get the appropriate /iua 
number equation N 



from the 



E 

k 



Indepen- 



dent optimization of /i^ar has a negligible effect on the 
energy. For our fully projected wave function of Eq. |6j 
we avoid attaching physical meaning to the parameters 
Arf and ^^ar- They are taken simply as two indepen- 
dent degrees of freedom in the variational wave function. 
In order to characterize the quantum phases, the appro- 
priate order parameters that characterize the SF and AF 
long range order are calculated explicitly in the optimized 
wave functions. It can be seen that for the wave function 
of Ref. [IB] , when long range hopping parameters t' and 
t" of ^k are set to zero, the ak parametrization of the 
wave function becomes equivalent to that of Ref. [13, 
and Eq. [Bl Unlike the case of high Tc materials where 
the t' and t" terms are important to describe the correct 
band structure of the materials, in the context of cold 
atoms in optical lattices, these extra parameters can be 
set to zero. 



The quantum phase at half filling (zero hole doping) 
has long range antifcrromagnctic (AF) order [l^. For 
large enough system sizes, a broken SU(2) symmetric 
ground state becomes likely. In a ground state that 
has antifcrromagnctic long range order, we note that 
even at zero doping, the singlet correlation accounts for 
large part of the energy contribution [soj . Thus, we 
keep the singlet correlation mechanism of Eq. [5] and ad- 
just the fluctuations around the Neel state with a Jas- 
trow factor that explicitly breaks the SU(2) symmetry: 
J ^e"'-s^-!. Here, 



ruaf 



latt 



ieLi 



(7) 



|*SF,) = JVGVN\SFd)MF 



is the staggered magnetization operator with = 
l/2(ni_+i — Ui^-i). Also, Lj is the sub-lattice of index 
j. Aaf can be interpreted as a local magnetic field that 
breaks the SU(2) symmetry. Thus, for the zero doping 
case we are assuming three variational parameters: A^, 
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Aay, and fiyar- The improvement in energy due to this 
choice of the Jastrow factor is ^ 2% which is much larger 
than the typical statistical errors. It was found that the 
optimum value Aaf = 0.35 not only produces lower en- 
ergy but also the correct staggered magnetization [29j . 

2) Homogeneous solution with coexistence of 
antiferromagnetism and superfluidity: A homoge- 



neous state with both the d-wave superfluid and stag- 
gered magnetic long range order(denoted by \AF + SF^)) 
is based on the mean field solution [IB-ii [Hi [H . This 
state is favorable close to the zero hole doping a; = up 
to a critical x = Xc- While for Xc < x < x'^, a pure d-wave 
superfluid (SF^) phase is favored. The initial mean field 
\AF + SFd)MF state is obtained by solving the Hamilto- 
nian 



(k+Q),2 



h.c.]} 



ko- 



k(T 



(8) 



r 



where the k vector is restricted to the reduced Brillouin with = ~2t{cos{kx) + cos{ky)) — fj^^ar and m7v(Q) = 
zone and Q = (tt, tt). The diagonalization produces two magnetic variational parameter. The corresponding 



spin density wave bands denoted by subindices a and 6. 
The dispersion relations are given as 



^k + a+Q + v/(^k-s^k-Q)2 4-m7v(Q)2 



Eb{if^) = (9) 



quasiparticle operators are 

C^k,a,fT = ttkCk.cr + Cr/^kCk-l-Q.cr 

rfk,b,o- = -o-/3kCk,o- + ttkCk+q.o-- (10) 
From these, the ground state ansatz can be written as 



\^AF+SF,) = JVgVn\AF + SFd) 



MF 



JPg 



E 



N/2 



-k,a,-l 



k,&, 



d' 

-l"-k,6.-l 



|0). 



(11) 



Here, = Ak/ Ea{k) + ^{Ea{k)y 



and Bk 



-Ak/ Sfc(k) + ^(£b(k))2 + Al) with d-wave symmet- 
ric parametrization Ak = Ad{cos(kx) — cos{ky))/2. Vn 
projects the wave function into the subspace of total 
number of particles N ^ Ni + N2 with equal popula- 
tions A^i = N2. For this variational wave function, the 
optimization parameters are A^, mAr(Q), and ^var- 



In the hole doped regime, the contribution from 
particle-particle correlations is small and the hole-hole 
correlations dominate [25|. We assume a Jastrow fac- 
tor J of the form J ~ g^/^X^i.i '''''j^'j where the hole 
density at i-th site is hi = hi^\hi^i and the correlation 
function is parametrized by Vij = vl\vi — Vj\. We find an 



improvement in energy of the order 0.1% with w < 0. 
Here, the distance between sites Vij — \vi — rj| must be 
between the closest images across the tilted boundaries 
as we later discuss. 



III. DISCUSSION OF RESULTS 



In analogy with the high materials, we consider JJ 
values in the Eq. |4] centered at 13.33 (in units of €). For 
the case of the simple single band model (Eq. [ij this \J 
value is solely controlled by the optical lattice depth Vo 
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while for the model of Eq. [2l (7 is controlled by Vo as 
well as by Feshbach resonance(magnetic field). For our 
discussion, we assume the parametrization of Eq. [5] with 
Uo = 13.33 and Sg/t varied. For the non-resonant case, 
this is just a way of getting different values of U, while 
for the resonant case Sg/t has physical meaning. 



The quantum expectation values are evaluated by 
Monte Carlo multi-dimensional integration. In order to 
avoid numerical instabilities along the k.^. = ky points, 
we take a patch of lattice sites with the tilted bound- 
ary conditions as our system frame 0, . The size of 
the system is Niatt ~ L'^ + 1 sites, where L is an odd 
number. In this way, we also avoid the frustration of 
the antiferromagnetic phase (one particle of spin -l-l(-l) 
should be surrounded by 4 particles of spin -l(-l-l)) at 
the half filling and low hole doping limits. The pairing 
wave function in real space is obtained by Fourier trans- 
forming (f>k = Uk/vk- This function is continuous when 
crossing the tilted boundaries. 



The injection of holes to the antiferromagnetic(AF) 
phase at half filling destabilizes the AF and beyond a crit- 
ical hole doping density Xc, the AF order is completely 
destroyed giving way to the homogeneous d-wave super- 
fluid (SFd) (Eq. (5]). There are two possible scenarios 
of AF f-> SFrf transitions considered here: 1) a first or- 
der phase transition between AF and SF^ phases as a 
function of doping. 2) a second order phase transition 
from a homogeneous AF + SFd phase at low doping den- 
sities to the SFd phase. The energetics of this scenario 
is calculated using the wave function of Eq. [TTJ In both 
cases, the pure SF^ phase exists in the hole doping regime 
Xc < X < x'^. For X > x'^^ the system becomes a normal 
Fermi fluid (NFF) phase with none of the long range an- 
tiferromagnetic and superfluid correlations. In all cases, 
the s-wave symmetric superfluidity is energetically disfa- 
vored. 
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FIG. 2: (Color online) Superfluid long range order $ (dia- 
monds, continuous line) and antiferromagnetic order param- 
eter AIaf (dashed line) at different hole densities and cou- 
pling U (plots a-c). These figures correspond to the scenario 
with phase separation. Here, we assumed the full Hamilto- 
nian of Eq. 3] and the interaction as parametrized hy U = 
Uo/{l + 5g/tf with Uo = 13.33 in units of t (Eq. [5]). In the 
phase separated region Q < x < Xc (horizontal dotted line), 
the order parameters are assumed to depend linearly on x, 
so that (^{x) = <b{xc)x/xc and Maf{x) = Maf{Q) [1 — x/xc]- 
For comparison, we also show Maf{x) for the homogeneous 
phase scenario (circles) at the interaction given by 5g/t = 
(plot a). We find that the order parameters for the phase 
separation scenario and the homogeneous phase scenario are 
very close to each other. U = 9.26 < Uo {Sg/t — 0.2, plot b) is 
more favorable for superfiuidity with a larger region of phase 
coexistence, while the opposite is true for U = 53.32 > Uo 
{Sg/t = —0.5, plot c). There is a distinguishable difference 
for the optimizing superfiuid parameter between Sg/t — 
and Sg/t — —0.5 (plot d). However, there is only a small dif- 
ference for Ad between Sg/t — and Sg/t — 0.2 cases; thus we 
use the same parametrization for both. The optimum value 
of the antiferromagnetic parameter Aaf was ~ 0.35 at x = 0. 
It has a rather weak dependence on Sg/t (that is, on U). 



A. Magnetic and Superfluid order parameters 



In order to characterize the quantum phases, we cal- 
culate the following long range orders; d-wave order $ = 

lim </.(j). Here, we defined q^{j) = , h/N {J^ sj^.B,) 

with Bj = l/2(4+i4+.__i - 4_icl+^,+i).^ And the 
long range order for the staggered magnetization Maf = 
{niaf) (Eq. [71). As shown in the Fig. [2J $ has a re- 
gion of favorable d-wave paring at finite hole densities. 
At zero hole density, we found that the broken SU(2) 
symmetry wave function produces the correct finite size 
behavior of Maf = 0.75(2) in agreement with previous 



variational calculations [26|, |29[ for the considered sys- 
tem sizes {Niatt — 50 to 226). These long range orders, 
however, cannot provide information on the nature of the 
phase transition (first or second order). In fact, they are 
shown to be in close agreement for the phase separation 
and the homogeneous phase cases (Fig. [5]). We find that 
for U > Uo where the exchange is suppressed, the su- 
perfluid pairing is also suppressed (see Fig. [2]), while in 
the U < Uo case the relative mobility t/U of the par- 
ticles is enhanced and it becomes favorable for d-wave 
superfluidity. 
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FIG. 3; (Color online) Energy density per hole ehoie{x) as a 
function of the hole density x is shown for the full Hamiltonian 
model. The interaction strengths corresponding to 5g/t = 
-0.5 {U/t = 53.32), (U/t = Uo/t = 13.33) and 0.2 {U/t = 
9.26) cases are considered. We observe a minimum of ehoie{x) 
at a; = > signaling phase separation in the region < 
X < Xc. Xc increases monotonically with decreasing U while 
still remaining in the large U regime. The results are obtained 
on systems of Niatt = 9^ + 1 = 82 lattice sites. The unit of 
energy is t. 





Mixed 




d-wave 




A A Homogeneous 




5g/t = -0.5 




5g/t = _ 




Sg/1 = 0.2 



FIG. 4: (Color online) The ground state energies of Hamilto- 
nian in Eq. |4]as a function of doping x for various interaction 
strengths parametrized hy Sg/t: Sg/t = —0.5 {U/t = 53.32), 
{U/t = Uo/t = 13.33) and 0.2 {U/t = 9.26). Within the 
region defined by a: £ {Q,Xc), we also show the energies cor- 
responding to the first order (continuous line) and the sec- 
ond order (triangle) phase transition scenarios. The energies 
are shown to be degenerate (within < 0.01% error bars) for 
the mixed and homogeneous cases. SF^j smoothly transitions 
to the normal Fermi fluid NFF aX x > x'^ {x'^ ~ 0.4) when 
Ad — O.The unit of energy is t. 



B. Phase separation vs Coexistence at low hole 
doping 



In order to characterize the inhomogeneous mixture 
phase, the thermodynamic considerations are as follow- 
ing: From thermodynamic constraints, the energy per 
site of the ground state £{x) = E{x)/Niatt has to be 
a convex function d^£{x) / dx^ > 0. The range of x S 
(xq^xi) where d£{x)/dx = constant (or d^£{x)/dx^ = 
0) imphes that the number of holes can be varied while 
the chemical potential /i = d£{x)/dx is kept the same. 
This is a signature of the first order phase transition. In 
this regime, the system is a mixture of phase / at Xf) and 
phase // at xi. In our case, xo = Q and phase / has 
antiferromagnetic order(AF) while Xc = xi and phase // 
is the d-wave superfluid(SF(i). In order to check whether 
there is an interval defined by a; = and Xc > 0, we 
should check for the flatness of £{x) [s^]- However, be- 
cause we lack the knowledge of the exact mixed phase 
ground state at all x we use the Maxwell construction: 
we calculate the energy per hole eh{x) and check for the 
existence of the minimum 1351 of 
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FIG. 5: (Color online) Condensation energy per site at differ- 
ent values of 5g/t for the Hamiltonian of Eq. |4] The unit of 
energy is t. Hole densities x between and Xc are considered. 
In this = corresponds to the pure antiferromagnet 

while X = Xc corresponds to the pure superfluid phase. For 
< a; < a;c, the condensation energy is that of the homoge- 
neous solution. 



efi{x) 



£iiix)-£iiO) 



(12) 



where £jjj{x) are the energy densities of the pure phases 
(AF and SF^^). This is equivalent to finding Xc such that 
^^dx^'' Uc = ^"'^^"^^^'^^^ . For all cases considered here, 
Ch {x) has a minimum (Fig. and Maxwell construction 
with non zero Xc is possible. In order to avoid pair break- 
ing effects, we introduce two holes at a time by removing 



one particle of each species. From Fig. [3J we find that 
the critical doping strengths are: for 6g/t — 0, Xc ^ 0.09, 
for 6g/t = 0.2, Xc « 0.12, and for Sg/t = -0.5, Xc ~ 0.05. 
Correspondingly, the long range orders within the region 
< .T < a;c are proportional to the areas of the pure 
quantum phases; that is, the AF phase with a; = and 
SFd phase with x = Xc (see the caption of Fig. [5]). 
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C. Ground State Energies 



The ground state energies for the full Hamiltonian are 
shown in the Fig. 2] We notice that the mixed phase 
and the homogeneous phase energies are closely degen- 
erate within the region of a; S (0,Xc) in all of the con- 
sidered cases. The kinetic energy has no contribution at 
zero doping. In this case, the only contribution comes 
from exchange. As holes are injected to the system, ex- 
change of particles becomes suppressed while the kinetic 
contribution increases. In the normal phase regime be- 
yond x'^, most of the contribution to the energy comes 
from the kinetic term. The composition of energies for 
the t-J model is shown in the Fig. |51 



From the Fig. [S] we see a monotonically decreasing be- 
havior of the condensation energy Econd which is the dif- 
ference in energies between the normal state and the state 
with quantum coherence. Since Econd is at its largest for 
the pure antifcrromagnetic phase, we can argue that once 
the system is phase separated it would not transition into 
the homogeneous phase. This would be true even when 
the superfluid portion of the system is destroyed by ther- 
mal fluctuations as long as the antifcrromagnetic portion 
retains its phase coherence. Thus, we argue that the 
phase separation scenario is more robust than the homo- 
geneous coexistence. 



in the two dimensional optical lattices. We have shown 
that the interaction strength U has a direct impact on 
the shape of the phase diagram. In the laboratory, the 
optical lattice is usually placed within a smoothly vary- 
ing harmonic trap. Then, the usual shell structure is ex- 
pected to appear. We have shown that the mixed phase 
is energetically allowed and possibly more robust than 
the homogeneous phase. 



There are several experimental techniques to measure 
the quantum correlations: for example, radio frequency 
techniques to directly measure the excitation spectrum 
of the Fermi gas are being developed (see Ref . [S^ ) • The 
measurement of the correlations in the noise could also 
be used, as it was applied to the detection of the Mott 
insulator phase (see Ref. [ill). Also, Bragg scattering 
is a scheme that allows the detection of the nodal points 
in the momentum space 



37| . A connection of the single 



band model to the resonant regime was made through 
perturbation theory, although a more elaborate multi- 
band model might be required for improved description. 
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D. Comparison of the Full vs t-J models 

V. APPENDIX 



In the Figs. [S] and [71 a comparison of the full and the t- 
/ Hamiltonians at different values of 5g/t is given. Both 
models converge at the zero doping x = since three- 
site hopping terms are density suppressed. However, at 
finite doping, differences arise: as seen in the Fig. [SJ the 
e-hoieix) tends to have a slightly larger curvature leading 
to a clear definition of Xc- This tendency is more no- 
ticeable for 5g/t > ([/ < Uo)- But, for U « Uo, the 
validity of the perturbative Hamiltonian (Eq. [3] and |3]) 
is being compromised. Here, the critical Xc seems to be 
identical for these two models. 



IV. CONCLUDING REMARKS 



We have investigated the appearance of different quan- 
tum phases in a system of strongly interacting fermions 



For the resonant fermions with kpa large, multi-band 
effects have to be included. A model by DuanfTol. [20j 
accounts them by an effective single-band Hubbard-type 
model. This approach physically corresponds to taking 
into account only the lowest energy but exact one-particle 
and two-particle states on each lattice site. For a single- 
particle state, the corresponding wave function is sim- 
ply the ground sate wave function in the local potential 
well. In the strongly interacting limit, the wave func- 
tion of a two-particle state includes contributions of all 
excited states and does not reduce to a properly sym- 
metrized product of the two ground state wave func- 
tions. As a result, the overlap of the wave functions for 
nearest neighbor sites and, therefore, the hopping am- 
plitude between them will depend on a particular den- 
sity configuration (see Fig. [1]): The hopping amplitude 
^01 = tio = t for the hop from a singly occupied site to an 
empty site, will be different from the hopping amplitudes 
til = to2 = t2o = t + Sg for a hop from a singly occu- 
pied site to another singly occupied one and for a reverse 
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parameter 8g/t = 



parameter 8g/t = 0.2 



parameter 5g/t = -0.5 



A 41-J 
o OFull 
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FIG. 6: (Color online) Comparison of the energy density per hole ehoie(x) for the full and the t-J Hamiltonians at different 
values of 5g/t: 5g/t = -0.5 {U /t = 53.32), {U /t = Uo/t = 13.33) and 0.2 {U/t = 9.26). For the full Hamiltonian, the 
minimum of ehoie(x) is more clearly defined signaling greater tendency for phase separation. The unit of energy is t. 



parameter 8g/t = 



parameter 5g/t = 0.2 



parameter 5g/t = -0.5 
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FIG. 7: (Color online) Comparison of the ground state energies for the full and the t-J Hamiltonians at different values of 5g/t: 
5g/t = -0.5 {U/t = 53.32), {U/t = Uo/t = 13.33) and 0.2 {U/t = 9.26). The full Hamiltonian includes three-site hopping 
terms (Eq. |4]) that lower the energy relative to the t-J model. At zero hole doping, three-site terms are suppressed and both 
models converge. The unit of energy is t. 



o - - -o Total 
A AKinetic 
*---*Exchanae 



act two-particle wave functions on nearest-neighbor sites 
and, hence, will be different from both toi and tn = to2- 
As a result, the Hamiltonian of the extended Hubbard 
model reads 

7^ = - ^ [t + 5g{n^,ij + nj^s) + 5t{ni^anj^s)]c\,^Cj^a 

<i.j>a 

-^Uo'^n.i^+ini^-i, (13) 



FIG. 8: (Color online) Doping dependence of the total, kinetic 
and exchange energies for the t-J model. The unit of energy 
is t. 



where the density assisted hopping terms in the kinetic 
energy is characterized by 6g and 5t. 



hop from a doubly occupied site to an empty one, re- 
spectively. This is because the amplitude toi contains 
the overlap between the two ground state wave func- 
tions on the nearest-neighbor sites, while the amplitudes 
ill = ^02 are determined by the overlap of the exact 
two-particle wave function on a site with the product of 
the ground state wave functions on the same and the 
nearest-neighbor sites. In turn, the hopping amplitude 
ti2 =^21 ~ t + 5t for a hop from a doubly occupied site 
to a singly occupied site contains the overlap of the ex- 



The coefficients 5g/t and 5t/t depend on the charac- 
teristics of the atomic species such as the width (in mag- 
netic field) of the resonance AS, the background scat- 
tering length and the difference of the atomic mag- 
netic moments between the closed and the open scatter- 
ing channels A/i. In particular, 5g/t can be written as 
6g/t «_(0.38(Vo/£^r)°-^'^ -l) + a (see the Fig. 1 of the 
Ref. [19[). a is a weakly dependent quantity of Vo/Er 
and proportional to AB\abg\Afi. Thus, for a given 
atomic species and a resonant magnetic field, dg/t can 
be tuned over a range of the optical lattice depths. 
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We define the kinetic operators JCm as 

<j,i>cr 

<i.j>(T 

<i,j>cr 

IC-1 = -{t + Sg) ^ hi^sc\,jC-j^anj^g (14) 

<j,i>cr 

change the number of doubly occupied sites by m 
and obeys commutation relation [V, Km] = mUolCm with 



Here, we have used the usual definition hi^^ = 1 — rii^^. 
The full Gutzwiller projector can be shown to drop the 
interaction and the 5t dependent terms. The effective 
Hamiltonian after Schrieffer- Wolff transformation is then 
the same as Eq. |4] with the definition of the effective 
U = C/o/(l + 5g/iY- In a recent work [s^ on the one 
dimensional attractive fermionic gas with population im- 
balance, it was found that the pairing order gets en- 
hanced with 5g < 0, while the spin-spin correlation is 
suppressed. Obviously, due to the difference in the di- 
mensionality (ID vs 2D) and the nature of the interac- 
tion (attractive vs repulsive), an intuitive connection is 
rather complicated. 
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